template<class T>
void swap(T &a, T &b) {
  T c;
  c=a;
  a=b;
  b=c;
}

mdp_real EofP(mdp_real kappa_s, 
	    mdp_real r_s, 
	    mdp_real zeta, 
	    mdp_real alpha_1, 
	    mdp_real alpha_2, 
	    mdp_real *p) {
  mdp_real kappa_t=kappa_s/zeta;
  mdp_real m0=0.5/kappa_t-(1.0+3.0*r_s*zeta);
  mdp_real p1=0;
  mdp_real p2=0;
  mdp_real y;
  int   i;
  mdp_real cs=-2.0*alpha_1;
  mdp_real cp=0; // +2.0*alpha_2/r_s;

  printf("m0=%f\n", m0);

  if(r_s!=0.0) 
    cp=+8.0*alpha_2/r_s; // hope this is right..?
  else 
    cp=0;

  //  mdp_real m0=0.5*(1.0/kappa-8);
  
  p1=
    pow((1.0+cs)*sin(p[1])-cs*0.5*sin(2.0*p[1]),2)+
    pow((1.0+cs)*sin(p[2])-cs*0.5*sin(2.0*p[2]),2)+
    pow((1.0+cs)*sin(p[3])-cs*0.5*sin(2.0*p[3]),2);
  
  p2=
    2.0*(1.0+cp)*(3.0-cos(p[1])-cos(p[2])-cos(p[3]))-
    0.5*cp*(3.0-cos(2.0*p[1])-cos(2.0*p[2])-cos(2.0*p[3]));
  
  y=1.0+(pow(m0+0.5*r_s*zeta*p2,2)+zeta*zeta*p1)/(2.0*(1.0+m0+0.5*r_s*zeta*p2));
  return log(y+sqrt(y*y-1));
}

void show_pion(fermi_field &psi) {
  site               x(mylattice);
  mdp_real              precision, c2[box[0]];
  mdp_real p[4];
  mdp_complex c[TMAX], tmp;
  int     n[4];
  float   px[TMAX];
  Measure py[TMAX];
  Measure pa[2];
  Measure E, Z;
  float   c3x[TMAX*TMAX];
  Measure c3y[TMAX*TMAX];
  mdp_real   c3t[TMAX*TMAX];

  if(strcmp(prompt("<stdin>", "PION", "NO"),"YES")==0) {
    fermi_propagator S(mylattice, 3);
    S.param=psi.param;
    S.precision=precision;
    // set_cold(U);
    generate(S,U);
    float   cx[TMAX];
    Measure pion[TMAX];
    for(t=0; t<box[0]; t++) {
      cx[t]=t;
      pion[t].mean=pion[t].error=0;
    }
    forallsites(x) {
      for(a=0; a<4; a++)
	for(b=0; b<4; b++)
	  pion[x(0)].mean+=real(trace(S(x,a,b)*hermitian(S(x,a,b))));
    }
    printf("pion propagator\n");
    for(t=0; t<box[0]/2; t++) 
      printf("%i,\t%e\n", t, pion[t].mean);
    for(t=0; t<box[0]; t++) pion[t]=log(pion[t]);
    linear_fit(cx, pion, 2, box[0]/2-2,pa);
    printf("log of pion propagator\n");
    for(t=0; t<box[0]/2; t++) 
      printf("%i,\t%e,\t%e\n", t, pion[t].mean, pa[0].mean*t+pa[1].mean);
    for(t=box[0]/2; t<box[0]; t++) 
      printf("%i,\t%e,\t%e\n", t, pion[t].mean, pa[0].mean*(box[0]-t)+pa[1].mean);
    printf("pa[0]=%e +- %e\n", pa[0].mean, pa[0].error);
    printf("pa[1]=%e +- %e\n", pa[1].mean, pa[1].error);
  }

  int source_spin=(int) val(prompt("<stdin>", "SOURCE_SPIN", "0"));
  int source_color=(int) val(prompt("<stdin>", "SOURCE_COLOR", "0"));
  int sink_spin=source_spin;
  int sink_color=source_color;


  forallsites(x) 
    for(a=0; a<4; a++)
      for(i=0; i<3; i++) 
	psi(x,a,i)=0;
  x.set(x0[0],x0[1],x0[2],x0[3]); psi(x,source_spin,source_color)=1;

  if(strcmp(prompt("<stdin>", "RELOAD", "NO"),"NO")==0) {
    
    time=mpi.time();
    mul_invQ(chi,psi,U,1,precision);
    printf("time: %f\n", mpi.time()-time);
    chi.save("fermi_field.mdp");

    int px, py, pz;
    t=6; px=0; pz=0; a=0; i=0;
    for(py=0; py<box[2]; py++) {
      x.set(t,px,py,pz);
      printf("x=(%i, %i, %i, %i) psi(x,spin=%i,color=%i)=%e + i %e\n",
	     t,px, py, pz, a, i, real(chi(x,a,i)), imag(chi(x,a,i)));
    }
    for(t=0; t<box[0]; t++) c2[t]=0;
    forallsites(x) {
      t=x(0);
      c2[t]+=real(chi(x,sink_spin,sink_color));
    }
    
    printf("quark propagator\n");
    for(t=0; t<box[0]; t++) printf("%i,\t%e\n",t, c2[t]);
  
  }

  if(strcmp(prompt("<stdin>", "SPEEDOFLIGHT", "NO"),"YES")==0) {

    chi.load("fermi_field.mdp");

    int p_counter=0;
    for(n[1]=0; n[1]<box[1]/2; n[1]++)
      for(n[2]=0; n[2]<box[2]/2; n[2]++)
	for(n[3]=0; n[3]<box[3]/2; n[3]++) {
	  
	  p[1]=(2.0*Pi*n[1])/(box[1]);
	  p[2]=(2.0*Pi*n[2])/(box[2]);
	  p[3]=(2.0*Pi*n[3])/(box[3]);
	  
	  if(p[1]*p[1]+p[2]*p[2]+p[3]*p[3]>=0) {
	    
	    printf("%i %i %i\n", n[1], n[2], n[3]);
	    
	    for(t=0; t<box[0]; t++) c[t]=0;
	    
	    forallsites(x) {
	      t=(x(0) + box[0] - x0[0]) % box[0];
	      
	      tmp=chi(x,sink_spin,sink_color)*
		(exp(-I*(p[1]*x(1)+p[2]*x(2)+p[3]*x(3))));
	      
	      c[t]+=tmp;
	    }
	    
	    mpi.add(c, TMAX);
	    
	    // and now fitting ////////////////////////////////
	    
	    if(ME==0) {	
	      for(t=0; t<box[0]; t++) {
		px[t]=t;
		py[t].mean=log(abs(c[t]));
		py[t].error=0;
	      }
	      
	      linear_fit(px, py, fit_min, fit_max, pa);
	      
	      E=-1.0*pa[0];
	      Z=exp(pa[1]);
	      
	      printf("-----------DATA--------------\n");
	      printf("P = ( %i, %i, %i )\n", n[1], n[2], n[3]);
	      printf("-----------------------------\n");
	      printf("t\tC[t]\t\tC[2] (fit)\n");
	      printf("-----------------------------\n");
	      for(t=0; t<box[0]; t++)
		printf("%i,\t%e,\t%e\n", t, py[t].mean,
		       pa[0].mean*t+pa[1].mean);
	      printf("-------- ANALYSIS -----------\n");
	      printf("P = ( %i, %i, %i )\n", n[1], n[2], n[3]);
	      printf("E (theory) = %e\n", 
		     EofP(psi.param.kappa,
			  psi.param.r_s,
			  psi.param.zeta,
			  ((FermilabActionParameters*) psi.param.data)->alpha[2],
			  ((FermilabActionParameters*) psi.param.data)->alpha[3],
			  p));
	      printf("E = %e +- %e\n", E.mean, E.error);
	      printf("Z = %e +- %e\n", Z.mean, Z.error);
	      printf("-----------------------------\n\n\n");
	      
	      c3x[p_counter]=pow(p[1],2)+pow(p[2],2)+pow(p[3],2);
	      c3y[p_counter]=E*E;	    
	      c3t[p_counter]=
		pow(EofP(psi.param.kappa,
			 psi.param.r_s,
			 psi.param.zeta,
			 ((FermilabActionParameters*) psi.param.data)->alpha[2],
			 ((FermilabActionParameters*) psi.param.data)->alpha[4],
			 p),2);
	      p_counter++;
	    }
	  }
	}
    
    if(ME==0) {
      // sort all in p^2
      
      
      for(i=1; i<p_counter; i++)
	for(j=0; j<i; j++)
	  if(c3x[i]<c3x[j]) {
	    // printf("swapping...\n");
	    swap(c3x[i], c3x[j]);
	    swap(c3y[i], c3y[j]);
	    swap(c3t[i], c3t[j]);
	  }
      
      // average redundant copies
      
      
      for(i=0; i<p_counter-1; i++) {
	for(j=i+1; (c3x[i]==c3x[j]) && (j<p_counter); j++);
	// printf("i=%i j=%i\n", i,j);
	if(j-i>1) {
	  
	  for(k=i+1; k<j; k++) c3y[i]=c3y[i]+c3y[k];
	  c3y[i]=c3y[i]/(j-i);
	
	  for(k=i+1; k<j; k++) c3t[i]=c3t[i]+c3t[k];
	  c3t[i]=c3t[i]/(j-i);
	  
	  p_counter-=(j-i-1);
	  for(k=i+1; k<p_counter; k++) {
	    // printf("moving...\n");
	    c3x[k]=c3x[k+j-i-1];
	    c3y[k]=c3y[k+j-i-1];
	    c3t[k]=c3t[k+j-i-1];
	  }
	}
      }
      

      for(k=1; k<p_counter; k++)
	c3y[k]=(c3y[k]-c3y[0])/c3x[k];
      c3y[0].mean=1;
      c3y[0].error=0;
      for(k=1; k<p_counter; k++)
	c3t[k]=(c3t[k]-c3t[0])/c3x[k];
      c3t[0]=1;
      
      printf("p^2\t\t(E^2-E0^2)/p^2\t(theory)\n");
      for(i=0; i<p_counter; i++)
	printf("%e,\t%e,\t%e\n",
	       c3x[i],
	       c3y[i].mean,
	       c3t[i]);
      
      linear_fit(c3x, c3y, 0, p_counter, pa);
      printf("pa[0]=%e +- %e\n", pa[0].mean, pa[0].error);
      printf("pa[1]=%e +- %e\n", pa[1].mean, pa[1].error);
    }
    
  }

  mpi.close_wormholes();
  return 0;
};
